---

title: “Informe PISA” author: “Laura Martínez González de Aledo” date: “11/7/2020” output: html_document

El siguiente estudio correspone a un analisis del Informe PISA del año 2006. El conjunto de datos se ha construido utilizando la puntuación media en Ciencias por país del Programa para la Evaluación Internacional de Estudiantes, junto con el GNI per cápita (paridad del poder adquisitivo, dólares de 2005), el índice educativo, el índice de salud y el índice de desarrollo humano de la ONU (HDI).

Variables:

• Country
• Overall
• Issues
• Explain
• Evidence
• Interest
• Support
• Income
• Health
• Edu
• HDI

El objetivo es modelizar la relación entre la puntuación media (OSS) y el resto de variables, utilizando modelos de splines y GAM. Se debe realizar CV cuando se pueda.

Librerías

library(tidyverse) 
library(broom) # modelos en df
library(flextable) # Tablas formateadas
library(mgcv) # para estimar gam 
library(reshape2) # melt
library(magrittr) # Pipe operators
library(janitor) # Clean names
library(skimr) # == Summarize
library(imputeTS) # para cambiar los valores nulos por la media
library(PerformanceAnalytics) # grafico de correlaciones
library(corrplot) # para ver outliers
library(gam) # para calcular el gam
library(rsample) # para el split
library(boot) # cv para modelos glm() 

Datos

En primer lugar importo la tabla, a continuación elimino las columnas no voy a utilizar. Y por último hago una limpieza de tabla, borrando los duplicados y cambiando los Na´s por la media en vez de eliminarlos porque perderíamos información.

pisa2006 <- read.csv("pisasci2006.csv")

pisa2006 %<>% clean_names()
colnames(pisa2006)  # Ahora los nombres de las columnas esta en minuscula

##  [1] "country"  "overall"  "issues"   "explain"  "evidence" "interest"
##  [7] "support"  "income"   "health"   "edu"      "hdi"

pisa2006 %<>% distinct(country,.keep_all= TRUE) # duplicados

summarise_all(pisa2006, funs(sum(is.na(.)))) # contamos valores nulos

##   country overall issues explain evidence interest support income health edu
## 1       0       8      8       8        8        8       8      4      4   6
##   hdi
## 1   6

pisa2006 <- na_mean(pisa2006) # cambiamos los NaN por la media

datos_pisa <- select(pisa2006, -issues, -explain, -evidence)
# quitamos las variables categoricas

head(datos_pisa) # vemos que nuestra no varia porque no teniamos duplicados

##      country  overall interest  support income health      edu       hdi
## 1    Albania 473.1404 528.1579 512.1754  0.599  0.886 0.716000 0.7240000
## 2  Argentina 391.0000 567.0000 506.0000  0.678  0.868 0.786000 0.7730000
## 3  Australia 527.0000 465.0000 487.0000  0.826  0.965 0.978000 0.9200000
## 4    Austria 511.0000 507.0000 515.0000  0.835  0.944 0.824000 0.8660000
## 5 Azerbaijan 382.0000 612.0000 542.0000  0.566  0.780 0.799322 0.8060169
## 6    Belgium 510.0000 503.0000 492.0000  0.831  0.935 0.868000 0.8770000

EDA

# libreria skimr
skim(datos_pisa)

Data summary

Name	datos_pisa
Number of rows	65
Number of columns	8
_______________________
Column type frequency:
character	1
numeric	7
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
country	0	1	4	24	0	65	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
overall	0	1	473.14	51.06	322.00	436.00	487.00	511.00	563.00	▁▂▂▇▃
interest	0	1	528.16	46.62	448.00	503.00	528.16	549.00	644.00	▆▇▇▃▂
support	0	1	512.18	24.40	447.00	497.00	512.18	527.00	569.00	▂▃▇▅▂
income	0	1	0.74	0.11	0.41	0.66	0.74	0.83	0.94	▁▂▆▇▅
health	0	1	0.89	0.06	0.72	0.84	0.89	0.94	0.99	▂▂▇▇▇
edu	0	1	0.80	0.11	0.54	0.74	0.80	0.87	0.99	▂▃▇▇▅
hdi	0	1	0.81	0.08	0.58	0.75	0.81	0.87	0.94	▁▃▆▇▆

Analisis de Correlaciones

Observamos las correlaciones para ver si hay problemas de multicolinealidad, dependiendo de si su correlación es lineal o no. Quito la variable Country porque no es una variable númerica

# librería corrplot
corrplot(cor(datos_pisa%>%
               select_at(vars(-country)),
             use = "complete.obs"),
         method = "circle", type = "upper")

# libreria PerformanceAnalytics
chart.Correlation(datos_pisa %>% 
                    select_at(vars(-country)),
                  histogram=TRUE, pch=12)

# Grafico interest
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = interest)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

# Grafico Support
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = support)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

# Grafico income
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = income)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

# Grafico health
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = health)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

# Grafico edu
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = edu)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

# Grafico hdi
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = hdi)) +
  layer(geom = "point",stat = "identity",position = "identity") +
  theme_bw() + theme(legend.key = element_blank())
baseplot1

Grados de Libertad

A continuación calculamos los grados de libertad de cada variable junto con el Cross Validation. Unicamente los calculo para las variables númericas.

datos_pisa <- select(pisa2006, -country, -issues, -explain, -evidence)

La función smooth.spline() de R permite ajustar smooth splines de forma sencilla, con la ventaja añadida de que el valor óptimo de smothness (λ) puede identificarse por cross-validation.

Overall Vs. Interest

plot(datos_pisa$interest, datos_pisa$overall, col='gray')
fit1 <- smooth.spline(datos_pisa$interest, datos_pisa$overall, df=16)
fit2 <- smooth.spline(datos_pisa$interest, datos_pisa$overall, cv=TRUE)
lines(fit1, col='red', lwd=2)
lines(fit2, col='blue', lwd=1)

fit2$df # 4.75

## [1] 4.750171

Overall Vs. Support

plot(datos_pisa$support, datos_pisa$overall, col='gray')
fit3 <- smooth.spline(datos_pisa$support, datos_pisa$overall, df=16)
fit4 <- smooth.spline(datos_pisa$support, datos_pisa$overall, cv=TRUE)
lines(fit3, col='red', lwd=2)
lines(fit4, col='blue', lwd=1)

fit4$df # 2.001

## [1] 2.001243

Overall Vs. Income

plot(datos_pisa$income, datos_pisa$overall, col='gray')
fit5 <- smooth.spline(datos_pisa$income, datos_pisa$overall, df=16)
fit6 <- smooth.spline(datos_pisa$income, datos_pisa$overall, cv=TRUE)
lines(fit5, col='red', lwd=2)
lines(fit6, col='blue', lwd=1)

fit6$df # 4.24

## [1] 4.244952

Overall Vs. Health

plot(datos_pisa$health, datos_pisa$overall, col='gray')
fit7 <- smooth.spline(datos_pisa$health, datos_pisa$overall, df=16)
fit8 <- smooth.spline(datos_pisa$health, datos_pisa$overall, cv=TRUE)
lines(fit7, col='red', lwd=2)
lines(fit8, col='blue', lwd=1)

fit8$df # 2.002

## [1] 2.002844

Overall Vs. Education

plot(datos_pisa$edu, datos_pisa$overall, col='gray')
fit9 <- smooth.spline(datos_pisa$edu, datos_pisa$overall, df=16)
fit10 <- smooth.spline(datos_pisa$edu, datos_pisa$overall, cv=TRUE)
lines(fit9, col='red', lwd=2)
lines(fit10, col='blue', lwd=1)

fit10$df # 2.00

## [1] 2.002385

Overall Vs. HDI

plot(datos_pisa$hdi, datos_pisa$overall, col='gray')
fit11 <- smooth.spline(datos_pisa$hdi, datos_pisa$overall, df=16)
fit12 <- smooth.spline(datos_pisa$hdi, datos_pisa$overall, cv=TRUE)
lines(fit11, col='red', lwd=2)
lines(fit12, col='blue', lwd=1)

fit12$df # 8.603

## [1] 8.603228

# suavizado de splines
fit2$lambda

## [1] 0.00465548

fit4$lambda

## [1] 81.82485

fit6$lambda

## [1] 0.007177517

fit8$lambda

## [1] 44.19162

fit10$lambda

## [1] 52.19016

fit12$lambda

## [1] 0.0002060338

Modelo GAM

# library(gam)
gam1 <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) +  s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603), data=datos_pisa)

summary(gam1)

## 
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) + 
##     s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi, 
##     8.603), data = datos_pisa)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.486  -8.202   3.138  11.680  51.442 
## 
## (Dispersion Parameter for gaussian family taken to be 596.3)
## 
##     Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 24091.7 on 40.402 degrees of freedom
## AIC: 620.1483 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## s(interest, 4.75)  1.000  90398   90398 151.5983 2.956e-15 ***
## s(support, 2.001)  1.000   2921    2921   4.8977   0.03260 *  
## s(income, 4.24)    1.000   4151    4151   6.9614   0.01177 *  
## s(health, 2.002)   1.000    745     745   1.2494   0.27027    
## s(edu, 2.002)      1.000    342     342   0.5731   0.45340    
## s(hdi, 8.603)      1.000   3320    3320   5.5670   0.02322 *  
## Residuals         40.402  24092     596                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                   Npar Df Npar F    Pr(F)   
## (Intercept)                                 
## s(interest, 4.75)     3.8 3.3660 0.020073 * 
## s(support, 2.001)     1.0 1.6190 0.210521   
## s(income, 4.24)       3.2 6.2341 0.001094 **
## s(health, 2.002)      1.0 1.4042 0.243012   
## s(edu, 2.002)         1.0 1.7007 0.199599   
## s(hdi, 8.603)         7.6 1.2307 0.307436   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(gam1, se=TRUE, col='blue')

gam2 <- gam(overall~ s(interest, 4.75) + s(support, 2.001)+  s(health, 2.002) + s(hdi, 8.603), data=datos_pisa)

summary(gam2)

## 
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) + 
##     s(health, 2.002) + s(hdi, 8.603), data = datos_pisa)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -81.927 -11.914   3.562  15.920  48.927 
## 
## (Dispersion Parameter for gaussian family taken to be 720.8222)
## 
##     Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 33622.08 on 46.6441 degrees of freedom
## AIC: 629.3297 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## s(interest, 4.75)  1.000 103182  103182 143.1452 8.117e-16 ***
## s(support, 2.001)  1.000   2016    2016   2.7975   0.10111    
## s(health, 2.002)   1.000   4341    4341   6.0220   0.01792 *  
## s(hdi, 8.603)      1.000    417     417   0.5788   0.45061    
## Residuals         46.644  33622     721                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                   Npar Df Npar F   Pr(F)  
## (Intercept)                               
## s(interest, 4.75)     3.8 3.7802 0.01094 *
## s(support, 2.001)     1.0 2.1137 0.15267  
## s(health, 2.002)      1.0 0.9677 0.33050  
## s(hdi, 8.603)         7.6 1.5333 0.17476  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(gam2, se=TRUE, col='blue')

gam3 <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(edu, 2.002), data=datos_pisa)

summary(gam3)

## 
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) + 
##     s(edu, 2.002), data = datos_pisa)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.389 -10.251   5.461  12.100  64.652 
## 
## (Dispersion Parameter for gaussian family taken to be 743.8547)
## 
##     Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 41095.59 on 55.2468 degrees of freedom
## AIC: 625.1709 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                       Df Sum Sq Mean Sq  F value  Pr(>F)    
## s(interest, 4.75)  1.000 103773  103773 139.5065 < 2e-16 ***
## s(support, 2.001)  1.000   2632    2632   3.5383 0.06524 .  
## s(edu, 2.002)      1.000   2821    2821   3.7930 0.05656 .  
## Residuals         55.247  41096     744                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                   Npar Df Npar F    Pr(F)   
## (Intercept)                                 
## s(interest, 4.75)     3.8 4.4461 0.004134 **
## s(support, 2.001)     1.0 2.6246 0.110885   
## s(edu, 2.002)         1.0 1.4429 0.234867   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(gam3, se=TRUE, col='blue')

ANOVA

Realizamos el test anova para comparar los tres modelos que hemos propuesto anteriormente.

anova.gam(gam1, gam2, gam3, test='F')

## Analysis of Deviance Table
## 
## Model 1: overall ~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) + 
##     s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603)
## Model 2: overall ~ s(interest, 4.75) + s(support, 2.001) + s(health, 2.002) + 
##     s(hdi, 8.603)
## Model 3: overall ~ s(interest, 4.75) + s(support, 2.001) + s(edu, 2.002)
##   Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1        65      24092                     
## 2        65      33622  0  -9530.4         
## 3        65      41096  0  -7473.5

Podemos comprobar que el que menor residuo tiene es el modelo 1 por lo que va a ser el modelo con el que vamos a trabajar.

Cross Validation

Una vez escogido el modelo, vamos a dividir nuestra base de datos entre muestra para test y muestra para entrenamiento, vamos a proceder a introducir nuestro modelo en el test para saber como predice.

Validación Simple:

Dividimos aleatoriamente las observaciones disponibles en dos grupos, uno se emplea para entrenar al modelo y otro para evaluarlo.

# library(rsample)
set.seed(123)
pisa_split <- initial_split(datos_pisa, prop =.7, strata = "overall")
pisa_train <- training(pisa_split)
pisa_test <- testing(pisa_split)

Tenemos la base de datos dividida en 70/30, y vamos a proceder a introducir nuestro modelo en el test para saber como predice.

gam_train <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) +  s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603), data=pisa_train)

plot(gam_train, se=TRUE, col='red')

summary(gam_train)

## 
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) + 
##     s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi, 
##     8.603), data = pisa_train)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.8039  -8.0293   0.5306   8.2329  27.3118 
## 
## (Dispersion Parameter for gaussian family taken to be 332.6423)
## 
##     Null Deviance: 108148.2 on 46 degrees of freedom
## Residual Deviance: 7451.927 on 22.4022 degrees of freedom
## AIC: 422.6816 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## s(interest, 4.75)  1.000  59311   59311 178.3015 3.793e-12 ***
## s(support, 2.001)  1.000   1418    1418   4.2632   0.05071 .  
## s(income, 4.24)    1.000   8218    8218  24.7065 5.387e-05 ***
## s(health, 2.002)   1.000   1921    1921   5.7749   0.02496 *  
## s(edu, 2.002)      1.000     47      47   0.1422   0.70966    
## s(hdi, 8.603)      1.000    287     287   0.8622   0.36302    
## Residuals         22.402   7452     333                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                   Npar Df Npar F    Pr(F)   
## (Intercept)                                 
## s(interest, 4.75)     3.8 3.1114 0.037830 * 
## s(support, 2.001)     1.0 0.9666 0.336116   
## s(income, 4.24)       3.2 6.1667 0.002732 **
## s(health, 2.002)      1.0 5.2212 0.032069 * 
## s(edu, 2.002)         1.0 0.1249 0.727573   
## s(hdi, 8.603)         7.6 2.5807 0.038446 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Leave One Out Cross-Validation (LOOCV)

La función cv.glm() calcula el error de predicción mediante cross-validation.

# library(boot)
cv <- cv.glm(data=datos_pisa, glmfit = gam1)
sqrt(cv$delta)

## [1] 35.19908 35.05841

K-fold Cross-Validation

Igual que el modelo anterior sólo que especificando K.

set.seed(123)
cv_error <- cv.glm(data=datos_pisa, glmfit = gam1, K=10)
cv_error$delta

## [1] 1250.154 1189.006

Cuando se especifica el número de grupos K a emplear en la validación, la función cv.glm() devuelve dos resultados, uno con corrección de continuidad y otro sin ella.

Referencias:

Smoothing Splines: https://rpubs.com/Joaquin_AR/250069